Table of Contents

  • Машинное обучение и предобработка данных
    • Постановка задачи.
      • Ход работы.
    • Исследование.
      • Загрузка, импорты. Описание данных. EDA.
      • Дубликаты и пропуски.
      • Анализ признаков, отбор, сокращение.
        • Корреляционный анализ.
        • Проверка стационарности.
        • Анализ информативности
        • Критерий Chi2
        • PCA
        • Графический анализ.
    • Обучение моделей.
      • Константная модель.
      • Логистическая регрессия.
      • K-ближайших соседей.
      • Случайный лес
    • Результаты на тестовой выборке, выводы, альтернативы.
      • Значимость признаков и модификация лучшей модели.
      • Матрица ошибок лучшей модели.
      • Некоторые альтернативы и компромиссы.
        • Бустинг.
        • Доп. логистическая регрессия.
        • Выбор порога threshold полученной модели.

Машинное обучение и предобработка данных
¶

выполнил: Морозов Е.А.

Постановка задачи.¶

В данной работе на примере одной из задач машинного обучения рассматривались различные способы: предобработки данных, отбора и анализа признаков, уменьшения размерности признакового пространства, подготовки данных перед их подачей в модель и популярные алгоритмы машинного обучения. Подход не претендует на полноту, он носит демонстративный характер для иллюстрации определённой стадии развития проф. навыков автора.

Источник данных (kaggle loan defaulter dataset): https://www.kaggle.com/datasets/gauravduttakiit/loan-defaulter/application_data.csv

Ход работы.¶

  1. Загрузка данных.

  2. Предобработка данных.

  3. Снижение размерности для визуализации данных.

  4. Обучение моделей классификации. Вывод и анализ метрик.

  5. Выводы.

Исследование.¶

Загрузим и проанализируем данные.

Датасет на странице соревнования представлен в разделённом 2 части виде, использовались лишь файлы "application_data.csv" и "columns_description.csv". Данные представлены в табличном виде, описание и анализ содержащихся признаков представлены ниже.

Для обучения моделей количество признаков избыточно и потребует большого объёма предообработки (от заполнения пропусков и исправления ошибок до кодирования переменных и анализа значимости), что растянет время подготовки и сделает этап отбора сложным в смысле написания большого конвейера предобработки. Нужно выработать подходы к отбору признаков для сокращения их чилса.

Загрузка, импорты. Описание данных. EDA.¶

Сделаем импорты, установим компоненты для работы, опишем и проанализируем наблюдаемые данные.

In [1]:
# установка недостающих библиотек и расширений
#!pip install ipywidgets
#!jupyter nbextension enable --py widgetsnbextension
#!pip install CatBoost
In [2]:
# импортируем библиотеки:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
sns.set_style('whitegrid')
np.set_printoptions(precision=4)
from random import sample

from sklearn.compose import ColumnTransformer, make_column_selector
from sklearn.decomposition import PCA
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import VarianceThreshold, mutual_info_classif, SelectKBest, chi2

from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.inspection import permutation_importance
from sklearn.metrics import f1_score, roc_auc_score, roc_curve, precision_score, recall_score,\
                            precision_recall_curve, accuracy_score, auc, classification_report, \
                            confusion_matrix, ConfusionMatrixDisplay

from sklearn.model_selection import train_test_split, RandomizedSearchCV, cross_val_score, GridSearchCV
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder

# отключим предупреждения:
import warnings
warnings.simplefilter('ignore')
pd.options.mode.chained_assignment = None # чтобы избавиться от конфликта pandas и sklearn

# зададим состояние генератора случайных чисел для исползуемых объектов:
state = np.random.RandomState(606)
#np.random.seed(606) - не работал как хотелось, не понял почему

Таблица изначально объемная с большим количеством объектов и признаков, в этой связи описание данных приводится не в 1-ой части работы, а из прилагающегося к датасету файла с описаниями колонок таблицы.

In [3]:
pd.read_csv('./4_2_data/columns_description.csv')
Out[3]:
Unnamed: 0 Table Row Description Special
0 1 application_data SK_ID_CURR ID of loan in our sample NaN
1 2 application_data TARGET Target variable (1 - client with payment diffi... NaN
2 5 application_data NAME_CONTRACT_TYPE Identification if loan is cash or revolving NaN
3 6 application_data CODE_GENDER Gender of the client NaN
4 7 application_data FLAG_OWN_CAR Flag if the client owns a car NaN
... ... ... ... ... ...
155 209 previous_application.csv DAYS_FIRST_DUE Relative to application date of current applic... time only relative to the application
156 210 previous_application.csv DAYS_LAST_DUE_1ST_VERSION Relative to application date of current applic... time only relative to the application
157 211 previous_application.csv DAYS_LAST_DUE Relative to application date of current applic... time only relative to the application
158 212 previous_application.csv DAYS_TERMINATION Relative to application date of current applic... time only relative to the application
159 213 previous_application.csv NFLAG_INSURED_ON_APPROVAL Did the client requested insurance during the ... NaN

160 rows × 5 columns

In [4]:
# загрузим датафрейм
df = pd.read_csv('./4_2_data/application_data.csv')
print(df.info())
df.sample(5)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 307511 entries, 0 to 307510
Columns: 122 entries, SK_ID_CURR to AMT_REQ_CREDIT_BUREAU_YEAR
dtypes: float64(65), int64(41), object(16)
memory usage: 286.2+ MB
None
Out[4]:
SK_ID_CURR TARGET NAME_CONTRACT_TYPE CODE_GENDER FLAG_OWN_CAR FLAG_OWN_REALTY CNT_CHILDREN AMT_INCOME_TOTAL AMT_CREDIT AMT_ANNUITY ... FLAG_DOCUMENT_18 FLAG_DOCUMENT_19 FLAG_DOCUMENT_20 FLAG_DOCUMENT_21 AMT_REQ_CREDIT_BUREAU_HOUR AMT_REQ_CREDIT_BUREAU_DAY AMT_REQ_CREDIT_BUREAU_WEEK AMT_REQ_CREDIT_BUREAU_MON AMT_REQ_CREDIT_BUREAU_QRT AMT_REQ_CREDIT_BUREAU_YEAR
282580 427314 0 Cash loans M N Y 0 202500.0 1078200.0 31653.0 ... 0 0 0 0 0.0 0.0 0.0 0.0 1.0 8.0
209933 343287 0 Cash loans M Y N 0 58500.0 95940.0 10822.5 ... 0 0 0 0 NaN NaN NaN NaN NaN NaN
51825 160018 0 Cash loans F N Y 0 112500.0 472500.0 46863.0 ... 0 0 0 0 0.0 0.0 0.0 0.0 0.0 2.0
77103 189406 0 Cash loans F N N 1 135000.0 381528.0 14382.0 ... 0 0 0 0 0.0 0.0 0.0 0.0 0.0 1.0
121059 240356 0 Cash loans F N Y 0 90000.0 239850.0 23494.5 ... 0 0 0 0 0.0 0.0 0.0 0.0 1.0 2.0

5 rows × 122 columns

Дубликаты и пропуски.¶

Таблица с данными весьма объёмная: 307.5К записей в >120 столбцах - содерижит как числовые значения, так и категориальные. Содержание пропусков и дубликатов рассматривается ниже и выполнено прежде остальной обработки, поскольку с отсутствующими данными никакие модели машинного обучения адекватно работать не в состоянии. Количество признаков явно избыточно для обучения: проявится тяга к переобучению, случайному характеру предсказаний на тестовой выборке и росту срока обучения из-за "проклятия размерности".

In [5]:
# проверим количество дубликатов
df.duplicated().sum()
Out[5]:
0
In [6]:
# проверим пропуски
nan_df = pd.DataFrame(df.isna().sum().sort_values(ascending=False), columns=['NaN_count'])
nan_df['ratio_%'] = 100*np.round(nan_df['NaN_count']/df.shape[0], 3)
nan_df.index.set_names('columns', inplace=True)

# примем процент пропусков, который будем считать критичным
top_nan_df = nan_df[nan_df['ratio_%']>10]
print(f"Столбцы таблицы данных с содержанием пропусков > 10% ({top_nan_df.shape[0]} шт.):")
top_nan_df

Столбцы таблицы данных с содержанием пропусков > 10% (57 шт.):
Out[6]:
NaN_count ratio_%
columns
COMMONAREA_MEDI 214865 69.9
COMMONAREA_AVG 214865 69.9
COMMONAREA_MODE 214865 69.9
NONLIVINGAPARTMENTS_MODE 213514 69.4
NONLIVINGAPARTMENTS_AVG 213514 69.4
NONLIVINGAPARTMENTS_MEDI 213514 69.4
FONDKAPREMONT_MODE 210295 68.4
LIVINGAPARTMENTS_MODE 210199 68.4
LIVINGAPARTMENTS_AVG 210199 68.4
LIVINGAPARTMENTS_MEDI 210199 68.4
FLOORSMIN_AVG 208642 67.8
FLOORSMIN_MODE 208642 67.8
FLOORSMIN_MEDI 208642 67.8
YEARS_BUILD_MEDI 204488 66.5
YEARS_BUILD_MODE 204488 66.5
YEARS_BUILD_AVG 204488 66.5
OWN_CAR_AGE 202929 66.0
LANDAREA_MEDI 182590 59.4
LANDAREA_MODE 182590 59.4
LANDAREA_AVG 182590 59.4
BASEMENTAREA_MEDI 179943 58.5
BASEMENTAREA_AVG 179943 58.5
BASEMENTAREA_MODE 179943 58.5
EXT_SOURCE_1 173378 56.4
NONLIVINGAREA_MODE 169682 55.2
NONLIVINGAREA_AVG 169682 55.2
NONLIVINGAREA_MEDI 169682 55.2
ELEVATORS_MEDI 163891 53.3
ELEVATORS_AVG 163891 53.3
ELEVATORS_MODE 163891 53.3
WALLSMATERIAL_MODE 156341 50.8
APARTMENTS_MEDI 156061 50.7
APARTMENTS_AVG 156061 50.7
APARTMENTS_MODE 156061 50.7
ENTRANCES_MEDI 154828 50.3
ENTRANCES_AVG 154828 50.3
ENTRANCES_MODE 154828 50.3
LIVINGAREA_AVG 154350 50.2
LIVINGAREA_MODE 154350 50.2
LIVINGAREA_MEDI 154350 50.2
HOUSETYPE_MODE 154297 50.2
FLOORSMAX_MODE 153020 49.8
FLOORSMAX_MEDI 153020 49.8
FLOORSMAX_AVG 153020 49.8
YEARS_BEGINEXPLUATATION_MODE 150007 48.8
YEARS_BEGINEXPLUATATION_MEDI 150007 48.8
YEARS_BEGINEXPLUATATION_AVG 150007 48.8
TOTALAREA_MODE 148431 48.3
EMERGENCYSTATE_MODE 145755 47.4
OCCUPATION_TYPE 96391 31.3
EXT_SOURCE_3 60965 19.8
AMT_REQ_CREDIT_BUREAU_HOUR 41519 13.5
AMT_REQ_CREDIT_BUREAU_DAY 41519 13.5
AMT_REQ_CREDIT_BUREAU_WEEK 41519 13.5
AMT_REQ_CREDIT_BUREAU_MON 41519 13.5
AMT_REQ_CREDIT_BUREAU_QRT 41519 13.5
AMT_REQ_CREDIT_BUREAU_YEAR 41519 13.5

Чтобы не заполнять пропуски и не разбираться в значениях излишне долго, столбцы с таким количеством пропусков удалим. Скорее всего это аггрегированные данные, полученные на основании информации из других столбцов таблицы, либо данные, полученные от слияния/объединения таблиц. Модели будут вполне способны приобрести обобщающую способность на оставшихся данных, поскольку их количество всё также велико. В случае необходимости можно будет извлечь признаки из столбца с датами отдельно.

In [7]:
# удаляем толбцы с большим кол-вом пропусков
df.drop(columns=top_nan_df.index, inplace=True)
df.isna().sum().sort_values(ascending=False).head(11)
Out[7]:
NAME_TYPE_SUITE             1292
OBS_30_CNT_SOCIAL_CIRCLE    1021
DEF_30_CNT_SOCIAL_CIRCLE    1021
OBS_60_CNT_SOCIAL_CIRCLE    1021
DEF_60_CNT_SOCIAL_CIRCLE    1021
EXT_SOURCE_2                 660
AMT_GOODS_PRICE              278
AMT_ANNUITY                   12
CNT_FAM_MEMBERS                2
DAYS_LAST_PHONE_CHANGE         1
FLAG_DOCUMENT_4                0
dtype: int64

Оставшиеся строки с пропусками можно как заполнить статистиками из соотв. распределений, либо вовсе удалить, поскольку даже если они и распределены по разным объектам, суммарное их количество не будет велико (порядка 5-6 тысяч записей из 305К).

In [8]:
# удаляем записи с пропусками
df.dropna(axis=0, inplace=True)
df.shape
Out[8]:
(304531, 65)

Анализ признаков, отбор, сокращение.¶

Выполнение задач этого пункта сопряжено с разведочным анализом данных, они будут пересекаться. Разделять выборук на признаки и целевой признак и на обучающую и тестовую части пока не станем, поскольку на данном этапе можно будет "отбросить" признаки из всей выборки сразу не опасаясь утечек, а после этого провести разделения.

Признаки могут демонстировать разной силы и природы связь между собой и с целевым признаком в зависимости от сочетаний из пары {категориальный тип данных; числовой тип данных}. Связь может быть как линейной, так и нелинейной, может отрицательной или остутствовать вовсе и в этой связи оценивать их полезность нужно отдельно и различными способами. Можно следовать согласно схеме, приведённой ниже.

features_selection.png

Используем некоторые пункты как основу для отбора признаков и проведём анализ:

  • корреляционный анализ числовых признаков,
  • дисперсионный анализ (не для целевого признака),
  • критерий chi2 для категориальных признаков и целевого признака,
  • анализ информативности,
  • компонентный анализ (как для иллюстрирования, так и для снижения кол-ва признаков).

Корреляционный анализ.¶

Быстрее всего рассмотреть линейные связи числовых признаков и удалить наиболее сильные чтобы избежать коллинеарности.

In [9]:
# рассмотрим матрицу корреляций (Pearson Correlation)
plt.figure(figsize=(12,10))
sns.heatmap(df.drop(columns=['SK_ID_CURR','TARGET']).corr()
            , cbar=True
            , annot=False
            , cmap=plt.cm.coolwarm_r)
plt.show()
In [10]:
corr_df = df.corr()

# функция для отбора признаков с высокой коррляцией

def correlation(dataset, threshold):
    col_corr = []
    corr_matrix = dataset.corr()
    for i in range(len(corr_matrix.columns)):
        for j in range(i):
            if corr_matrix.iloc[i, j] > threshold: # сравниваем не по модулю, т.к. отриц.
                                                   # корреляции могут нести инфомрацию для
                                                   # обучения моделей (предположительно)
                colname = corr_matrix.columns[i]   # получаем название столбца
                col_corr.append(colname)
    return col_corr

# выводим список сильно-скореллированных столбцов и удаляем их
f_correlated = correlation(corr_df, 0.75)
print('Признаки с наибольшими корреляциями:')
pd.Series(f_correlated, index=range(len(f_correlated)))

Признаки с наибольшими корреляциями:
Out[10]:
0                     AMT_ANNUITY
1                 AMT_GOODS_PRICE
2                 AMT_GOODS_PRICE
3                  FLAG_EMP_PHONE
4                 CNT_FAM_MEMBERS
5     REGION_RATING_CLIENT_W_CITY
6     LIVE_REGION_NOT_WORK_REGION
7         LIVE_CITY_NOT_WORK_CITY
8        OBS_60_CNT_SOCIAL_CIRCLE
9        DEF_60_CNT_SOCIAL_CIRCLE
10                FLAG_DOCUMENT_6
dtype: object
In [11]:
df.drop(columns=f_correlated, inplace=True)
df.shape
Out[11]:
(304531, 55)

Проверка стационарности.¶

После корр.анализа количество признаков вновь уменьшилось за счёт удаления столбцов с числовыми типами данных.

Однако также стоит проверить столбцы и на стационарность содержащихся значений, т.к. банальная логика подсказывает, что если в столбце напротив всех объектов содержится лишь одно неизменное значение или единицы раз встречается ещё несколько других значений, то сомнительно что подобный признак окажется полезным. Будем либо считать такие строки выбросными, либо удалять "стационарный" столбец целиком. Речь, конечно же, идёт не о целевом признаке.

In [12]:
# инициализация класса, для анализа zero-variance признаков
var_thres=VarianceThreshold(threshold=0)

# fit на числовые признаки и выводим список с True/False значениями
var_thres.fit(df.select_dtypes(include=['float','int']))
var_thres.get_support()   # столбцы с False стационарны
Out[12]:
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True, False,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True])
In [13]:
# отфильтровываем и удаляем столбцы
num_cols = df.select_dtypes(include=['float','int']).columns.values
constant_columns = [column for column in num_cols
                    if column not in df.select_dtypes(include=['float','int'])
                                       .columns[var_thres.get_support()]]
print(f'Стационарных признаков: {len(constant_columns)} шт. -', *constant_columns)
print('Признаков до удаления -', df.shape[1])
df.drop(columns=constant_columns, inplace=True)
print('Признаков после удаления -', df.shape[1])

Стационарных признаков: 1 шт. - FLAG_MOBIL
Признаков до удаления - 55
Признаков после удаления - 54

Анализ информативности¶

Полезность оставшихся признаков можно оценивать с помощью критерия *взаимной информации*

взаимная информация — статистическая функция двух случайных величин, описывающая количество информации, содержащейся в одной случайной величине относительно другой.

Определяется через энтропию и условную энтропию двух случайных величин как:

$ I(X;Y)=H(X)-H(X\mid Y)=H(X)+H(Y)-H(X,Y)$

Соотв. классы описаны в библиотеке sklearn, инициализируем их и обучим на столбцах датафрейма.

  • Однако для получения значений энтропии по всем столбцам текстовые значения категориальных переменных придётся перекодировать в числовые.

  • Преобразования выполняются пока что для всей выборки целиком, поскольку целью является выявление признаков, которые можно удалить как недостаточно информативные для обучения моделей.

  • Ограничиться использованием только прямого кодирования с пом. LabelEncoder() не вышло поскольку он не поддерживает 2D-массивы данных на входе и преназначен для работы с целевым признаком (1D-массив). Для обработки признаков использовался OrdinalEncoder().

In [14]:
# разделим признаки по категориям
cat_features = df.select_dtypes(include='object').columns.values     # категориальные
num_features = df.select_dtypes(include=['float','int']).columns.values    # числовые (только с плав.точкой, целые не масштабируем)
print(cat_features.shape, num_features.shape)

(11,) (43,)
In [15]:
# преобразуем признаки в ColumnTransformer
ct = ColumnTransformer([
     ("text_preprocess", OrdinalEncoder(), cat_features)
    #, ("num_preprocess", MinMaxScaler(), num_features)
], remainder='passthrough')

df_encoded = ct.fit_transform(df.drop(columns=['SK_ID_CURR', 'TARGET']))
df_encoded.shape
Out[15]:
(304531, 52)
In [16]:
# определим значения mutual information для перекодированных признаков относительно целевого
mutual_info = mutual_info_classif(
    pd.DataFrame(df_encoded, columns=df.columns[2:]), df['TARGET']
    , random_state=state, n_neighbors=30
)
In [17]:
# отсортируем полученные значения взаимной информации и отобразим на диаграмме
mutual_info = pd.Series(mutual_info, index=df.columns[2:])
mutual_info.sort_values(ascending=False).plot.bar(figsize=(20, 8))
plt.show()
  • На диаграмме наглядно показаны степени информативности признаков, поскольку малые величины взаимной информации легче воспринимаются в таком виде.
  • Для примера (и исходя из вида диаграммы) попробуем остановиться примерно на 30 наиболее информативных признаках с наибольшими значениями взаимной информации.
In [18]:
#top 30 признаков через SelectKBest, чтобы не выставлять и не определять пороги взятия вручную
sel_30_cols = SelectKBest(mutual_info_classif, k=30)
sel_30_cols.fit(df_encoded, df['TARGET'].values)
Out[18]:
SelectKBest(k=30,
            score_func=<function mutual_info_classif at 0x00000281FC43E940>)
In [19]:
#из атрибута .get_support() получаем метки является или нет признак значимым
sel_30_cols.get_support()
Out[19]:
array([False,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True, False,  True,  True, False, False,
        True,  True,  True,  True,  True,  True, False,  True, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False])
In [20]:
#по маске получаем наименования столбцов, учитвая что в df отбрасывались первые 2 столбца
top_30_cols = df[df.columns[2:]].columns[sel_30_cols.get_support()]
top_30_cols
Out[20]:
Index(['CODE_GENDER', 'FLAG_OWN_CAR', 'FLAG_OWN_REALTY', 'CNT_CHILDREN',
       'AMT_INCOME_TOTAL', 'AMT_CREDIT', 'NAME_TYPE_SUITE', 'NAME_INCOME_TYPE',
       'NAME_EDUCATION_TYPE', 'NAME_FAMILY_STATUS', 'NAME_HOUSING_TYPE',
       'REGION_POPULATION_RELATIVE', 'DAYS_BIRTH', 'DAYS_EMPLOYED',
       'DAYS_REGISTRATION', 'DAYS_ID_PUBLISH', 'FLAG_WORK_PHONE',
       'FLAG_CONT_MOBILE', 'FLAG_PHONE', 'FLAG_EMAIL', 'REGION_RATING_CLIENT',
       'HOUR_APPR_PROCESS_START', 'REG_REGION_NOT_LIVE_REGION',
       'REG_CITY_NOT_WORK_CITY', 'ORGANIZATION_TYPE', 'EXT_SOURCE_2',
       'OBS_30_CNT_SOCIAL_CIRCLE', 'DEF_30_CNT_SOCIAL_CIRCLE',
       'DAYS_LAST_PHONE_CHANGE', 'FLAG_DOCUMENT_3'],
      dtype='object')
In [21]:
#оставляем выбранные столбцы в качестве полезных признаков
df_reduced = pd.concat([df['TARGET'],df[top_30_cols]], axis=1)
df_reduced.head(3)
Out[21]:
TARGET CODE_GENDER FLAG_OWN_CAR FLAG_OWN_REALTY CNT_CHILDREN AMT_INCOME_TOTAL AMT_CREDIT NAME_TYPE_SUITE NAME_INCOME_TYPE NAME_EDUCATION_TYPE ... REGION_RATING_CLIENT HOUR_APPR_PROCESS_START REG_REGION_NOT_LIVE_REGION REG_CITY_NOT_WORK_CITY ORGANIZATION_TYPE EXT_SOURCE_2 OBS_30_CNT_SOCIAL_CIRCLE DEF_30_CNT_SOCIAL_CIRCLE DAYS_LAST_PHONE_CHANGE FLAG_DOCUMENT_3
0 1 M N Y 0 202500.0 406597.5 Unaccompanied Working Secondary / secondary special ... 2 10 0 0 Business Entity Type 3 0.262949 2.0 2.0 -1134.0 1
1 0 F N N 0 270000.0 1293502.5 Family State servant Higher education ... 1 11 0 0 School 0.622246 1.0 0.0 -828.0 1
2 0 M Y Y 0 67500.0 135000.0 Unaccompanied Working Secondary / secondary special ... 2 9 0 0 Government 0.555912 0.0 0.0 -815.0 0

3 rows × 31 columns

Корреляционный анализ и проверку стационарности повторно проводить не станем.

Критерий Chi2¶

Значимость категориальных признаков в задаче классификации можно проверить также с помощью критерия Fisher Score, проведя тест Chisquare. Эту оценку можно использовать для выбора признаков с наибольшими значениями тестовой статистики хи-квадрат, которая должна содержать только неотрицательные признаки, такие как логические значения или частоты относительно классов целевого признака.

Тест хи-квадрат измеряет зависимость между стохастическими переменными, поэтому использование этой функции «отсеивает» признаки, которые с наибольшей вероятностью не зависят от класса и, следовательно, не имеют значения для классификации. Статистика хи-квадрат обычно используется для проверки взаимосвязей между категориальными переменными. Она сравнивает наблюдаемое распределение различных классов целевых объектов Y среди различных категорий объектов с ожидаемым распределением целевых классов независимо от категорий объектов.

Поскольку перекодированные в числовые величины признаки у нас уже есть, то проведём также и данное тестирование:

  • критерий chi2 (Хи-квадрат) будет рассчитывать как значение статистики хи-квадрат, так и значение p-value.
  • будем отталкиваться от значений хи-квадрат – чем они больше, тем наиболее зависят признаки от таргета.
  • нулевая гипотеза в данном случае (провеяется по p-value): отсутствие связи между двумя переменными. При 0 соотв. будет показывать, что связь есть. Как и при значениях меньше порогового. Потому оставим столбцы, в которых p-value = 0 и меньше порогового.
In [22]:
# разделим признаки по категориям
cat_features_1 = df_reduced.select_dtypes(include='object').columns.values #категориальные
num_features_1 = df_reduced.select_dtypes(include=['float','int']).columns.values #числовые
print(cat_features_1.shape, num_features_1.shape)

(9,) (22,)
In [23]:
# преобразуем признаки в ColumnTransformer

ct_1 = ColumnTransformer([
     ("text_preprocess", OrdinalEncoder(), cat_features_1)
    #, ("num_preprocess", MinMaxScaler(), num_features_1)
], remainder='drop')

df_reduced_encoded = pd.DataFrame(
    ct_1.fit_transform(df_reduced.drop(columns=['TARGET']))
    , columns=cat_features_1
)
print(df_reduced_encoded.shape)

(304531, 9)
In [24]:
#выполним chi2 test
## chi2 возращает 2 значения: F-score и p-value
f_p_values=chi2(df_reduced_encoded, df_reduced['TARGET'])
f_p_values
Out[24]:
(array([6.1115e+02, 9.5505e+01, 4.2223e+00, 1.3310e+01, 9.2696e+02,
        4.7672e+02, 6.3285e-03, 2.5660e+02, 4.1864e+03]),
 array([6.2916e-135, 1.4750e-022, 3.9896e-002, 2.6393e-004, 1.3557e-203,
        1.1076e-105, 9.3659e-001, 9.4491e-058, 0.0000e+000]))
In [25]:
p_values = pd.Series(f_p_values[1], index=cat_features_1).sort_values(ascending=True)
p_values
Out[25]:
ORGANIZATION_TYPE       0.000000e+00
NAME_INCOME_TYPE       1.355740e-203
CODE_GENDER            6.291614e-135
NAME_EDUCATION_TYPE    1.107582e-105
NAME_HOUSING_TYPE       9.449137e-58
FLAG_OWN_CAR            1.474985e-22
NAME_TYPE_SUITE         2.639264e-04
FLAG_OWN_REALTY         3.989641e-02
NAME_FAMILY_STATUS      9.365937e-01
dtype: float64

  • Судя по мизерным величинам p-value у лидеров ("NAME_INCOME_TYPE", "CODE_GENDER", "NAME_EDUCATION_TYPE") они являются наиболее важными колонками по отношению к классам целевого признака.

  • При этом столбец "ORGANIZATION_TYPE" отмечен нулевой оценкой, т.е. противоречит 0-гипотезе, т.е. говорит о связи между признаками.

  • "NAME_FAMILY_STATUS" - обладает очень высоким значением p-value, большим чем любой общепринятый уровней значимости. Этот столбец, пожалуй, стоит удалить, в том числе и для того чтобы не считать данный этап бесполезным.

In [26]:
df_reduced.drop(columns=['NAME_FAMILY_STATUS'], inplace=True)
df_reduced.shape
Out[26]:
(304531, 30)

PCA¶

Сократив количество признаков за счёт прямого удаления из рассматриваемой выборки мы всё-таки продолжаем видеть большое их количество при столь же большом разнообразии, что на имеющемся объёме объектов делает задачу обучения через чур многомерной и наврядли позволит значимо превысить результат случайного классификатора.

  • Для сокращения размерности и увеличения плотности информации в признаках используется компонентный анализ, позволяющий перейти от одного признакового пространства в пространство меньшей размерности за счёт выделения новых признаков на основаниии распределений значений в имеющихся.

  • Особенно удобен он когда объекты нужно разделять линейно и при слабой интерпретируемости существующих признаков, для удаления же малозначимой информации при сокращении он также используется часто. Недостаток состоит лишь в том, что полученные признаки будут неитерпретируемыми, хотя и позволят проще проводить кластеризицию и разделения, обратные преобразования всю исходную информацию не восстановят.

  • По всей видимости стоит учесть также одновременное наличие числовых и категориальных величин в наборе данных. Самым логичным ходом видится их параллельная обработка с объединением в единый набор данных перед подачей в модель. Числовые признаки также стоит нормализовать перед проведение PCA, чтобы большие масштабы данных не угнетали малые в итоговом результате.

In [27]:
# набор данных пора разделить
X = df_reduced.drop(columns=['TARGET']).copy()
y = df_reduced['TARGET'].copy()

# рассмотрим балансы классов категориальных и дискретных числовых признаков

print(f'\nБаланс классов целевого признака.\n')

_y = pd.DataFrame(y.value_counts(normalize=True, ascending=False))
display(_y.style.format('{:.1%}'))

plt.figure(figsize=(5,5))
plt.pie(_y['TARGET'],labels=_y.index,autopct='%1.1f%%')
plt.show()

print(f'\nБаланс классов категориальных признаков c круговыми диаграммами.\n')

for col in df_reduced.columns:
    if (df_reduced[col].dtype == 'object'):
        print(f'\t- столбец {col}')
        _t = pd.DataFrame(df_reduced[col].value_counts(normalize=True, ascending=False))
        display(_t.style.format('{:.1%}'))

        plt.figure(figsize=(5,5))
        plt.pie(_t[col],labels=_t.index,autopct='%1.1f%%')
        plt.show()

Баланс классов целевого признака.
  TARGET
0 91.9%
1 8.1% 
Баланс классов категориальных признаков c круговыми диаграммами.

	- столбец CODE_GENDER
  CODE_GENDER
F 65.8%
M 34.2%
XNA 0.0% 
	- столбец FLAG_OWN_CAR
  FLAG_OWN_CAR
N 66.0%
Y 34.0% 
	- столбец FLAG_OWN_REALTY
  FLAG_OWN_REALTY
Y 69.5%
N 30.5% 
	- столбец NAME_TYPE_SUITE
  NAME_TYPE_SUITE
Unaccompanied 81.1%
Family 13.1%
Spouse, partner 3.7%
Children 1.1%
Other_B 0.6%
Other_A 0.3%
Group of people 0.1% 
	- столбец NAME_INCOME_TYPE
  NAME_INCOME_TYPE
Working 51.7%
Commercial associate 23.2%
Pensioner 18.0%
State servant 7.1%
Unemployed 0.0%
Student 0.0%
Businessman 0.0%
Maternity leave 0.0% 
	- столбец NAME_EDUCATION_TYPE
  NAME_EDUCATION_TYPE
Secondary / secondary special 71.1%
Higher education 24.2%
Incomplete higher 3.3%
Lower secondary 1.2%
Academic degree 0.1% 
	- столбец NAME_HOUSING_TYPE
  NAME_HOUSING_TYPE
House / apartment 88.7%
With parents 4.8%
Municipal apartment 3.6%
Rented apartment 1.6%
Office apartment 0.8%
Co-op apartment 0.4% 
	- столбец ORGANIZATION_TYPE
  ORGANIZATION_TYPE
Business Entity Type 3 22.1%
XNA 18.0%
Self-employed 12.5%
Other 5.4%
Medicine 3.6%
Business Entity Type 2 3.4%
Government 3.4%
School 2.9%
Trade: type 7 2.6%
Kindergarten 2.2%
Construction 2.2%
Business Entity Type 1 1.9%
Transport: type 4 1.8%
Trade: type 3 1.1%
Industry: type 9 1.1%
Industry: type 3 1.1%
Security 1.1%
Housing 1.0%
Industry: type 11 0.9%
Military 0.9%
Bank 0.8%
Agriculture 0.8%
Police 0.8%
Transport: type 2 0.7%
Postal 0.7%
Security Ministries 0.6%
Trade: type 2 0.6%
Restaurant 0.6%
Services 0.5%
University 0.4%
Industry: type 7 0.4%
Transport: type 3 0.4%
Industry: type 1 0.3%
Hotel 0.3%
Electricity 0.3%
Industry: type 4 0.3%
Trade: type 6 0.2%
Insurance 0.2%
Industry: type 5 0.2%
Telecom 0.2%
Emergency 0.2%
Industry: type 2 0.2%
Advertising 0.1%
Realtor 0.1%
Culture 0.1%
Industry: type 12 0.1%
Trade: type 1 0.1%
Mobile 0.1%
Legal Services 0.1%
Cleaning 0.1%
Transport: type 1 0.1%
Industry: type 6 0.0%
Industry: type 10 0.0%
Religion 0.0%
Industry: type 13 0.0%
Trade: type 4 0.0%
Trade: type 5 0.0%
Industry: type 8 0.0%
In [28]:
# классы несбалансированные, частично учтём это при разделении выборки
X_train, X_test, y_train, y_test = (
    train_test_split(X, y, train_size = 0.7
                     , random_state=state, shuffle=True
                     , stratify=pd.concat([y,
                                           X[['CODE_GENDER','FLAG_EMAIL','FLAG_OWN_CAR']]],
                                          axis=1))
)

X_train.shape, X_test.shape, y_train.shape, y_test.shape
Out[28]:
((213171, 29), (91360, 29), (213171,), (91360,))
In [29]:
# используя конвейеры подготовим данные для обучения и сожмём с пом. PCA

## описание для выбора категориальных и числовых данных
cat_features = make_column_selector(dtype_include=['object','int64'])
num_features = make_column_selector(dtype_exclude=['object','int64'])

## трансформеры для категориальных и числовых данных (параллельные)
ct_num = ColumnTransformer([
    ("num_preprocess", StandardScaler(), num_features)
], remainder='drop')

ct_cat = ColumnTransformer([
    ("cat_preprocess"
     , OrdinalEncoder(
         handle_unknown='use_encoded_value', unknown_value = -1
     )
     , cat_features)
], remainder='drop')

## преобразования данных из соотв. трансформмеров и сжатие PCA
num_pca_pipe = Pipeline([('num_prepr', ct_num), ('num_PCA', PCA(0.95))])
cat_pca_pipe = Pipeline([('cat_prepr', ct_cat),
                         #('cat_stand', MinMaxScaler()), #информ-ть масштабир. не повысило
                         ('cat_PCA', PCA(0.95))])

## объединение выходов с предыдущего шага в один набор данных для дальнейшего обучения
prepr_union = FeatureUnion([('numeric', num_pca_pipe),
                            ('categorical', cat_pca_pipe)])

## смаштабируем полученные признаки в ещё одном конвейре
pca_united_scaled = Pipeline([('features_union', prepr_union),
                              ('std_scaler', StandardScaler())])

print('''\nИспользуя конвейеры и трансформеры, подготовим данные для обучения моделей:
- параллельно обработаем числовые и категориальные признаки (нормализация, кодирование),
- сократим пространство признаков с помощью PCA
- объединим резульаты обработки категор. и числовых признаков в один набор данных

В качестве примера перехода от исходного набора данных к итоговому рассмотрим в сокращении:

\t1. исх.данные (после удаления столбцов с пропусками):''')
display(df.sample(3))

print('\n\t2. данные перед PCA (после удаления столбцов в рез-ате корр.анализа и теста Фишера):')
display(X.sample(3))
      
print('\n\t3. данные после PCA (0.95)')    
display(pd.DataFrame(pca_united_scaled.fit_transform(X_train)).sample(5))

Используя конвейеры и трансформеры, подготовим данные для обучения моделей:
- параллельно обработаем числовые и категориальные признаки (нормализация, кодирование),
- сократим пространство признаков с помощью PCA
- объединим резульаты обработки категор. и числовых признаков в один набор данных

В качестве примера перехода от исходного набора данных к итоговому рассмотрим в сокращении:

	1. исх.данные (после удаления столбцов с пропусками):
SK_ID_CURR TARGET NAME_CONTRACT_TYPE CODE_GENDER FLAG_OWN_CAR FLAG_OWN_REALTY CNT_CHILDREN AMT_INCOME_TOTAL AMT_CREDIT NAME_TYPE_SUITE ... FLAG_DOCUMENT_12 FLAG_DOCUMENT_13 FLAG_DOCUMENT_14 FLAG_DOCUMENT_15 FLAG_DOCUMENT_16 FLAG_DOCUMENT_17 FLAG_DOCUMENT_18 FLAG_DOCUMENT_19 FLAG_DOCUMENT_20 FLAG_DOCUMENT_21
37231 143126 0 Cash loans F N N 0 180000.0 485640.0 Unaccompanied ... 0 0 0 0 0 0 0 0 0 0
52646 160972 0 Cash loans F N N 0 157500.0 1293502.5 Unaccompanied ... 0 0 0 0 0 0 0 0 0 0
231955 368671 0 Cash loans F Y N 1 112500.0 1051245.0 Unaccompanied ... 0 0 0 0 0 0 0 0 0 0

3 rows × 54 columns

	2. данные перед PCA (после удаления столбцов в рез-ате корр.анализа и теста Фишера):
CODE_GENDER FLAG_OWN_CAR FLAG_OWN_REALTY CNT_CHILDREN AMT_INCOME_TOTAL AMT_CREDIT NAME_TYPE_SUITE NAME_INCOME_TYPE NAME_EDUCATION_TYPE NAME_HOUSING_TYPE ... REGION_RATING_CLIENT HOUR_APPR_PROCESS_START REG_REGION_NOT_LIVE_REGION REG_CITY_NOT_WORK_CITY ORGANIZATION_TYPE EXT_SOURCE_2 OBS_30_CNT_SOCIAL_CIRCLE DEF_30_CNT_SOCIAL_CIRCLE DAYS_LAST_PHONE_CHANGE FLAG_DOCUMENT_3
154402 F N Y 1 81000.0 450000.0 Unaccompanied State servant Secondary / secondary special House / apartment ... 2 9 0 0 Postal 0.282393 4.0 2.0 -214.0 1
274280 M Y Y 0 405000.0 585000.0 Unaccompanied Working Higher education House / apartment ... 2 14 0 0 Bank 0.641143 4.0 0.0 -1747.0 1
244999 M Y Y 0 202500.0 269550.0 Unaccompanied Commercial associate Secondary / secondary special House / apartment ... 2 9 0 1 Business Entity Type 2 0.181407 2.0 1.0 -406.0 1

3 rows × 29 columns

	3. данные после PCA (0.95)
0 1 2 3 4 5 6 7 8 9 10
39334 -0.240036 -0.662380 -0.074217 -0.755159 -0.899376 1.022755 1.164312 -0.313709 0.750766 0.194059 1.352621
91360 -0.039928 -0.356312 -0.140155 0.977006 -1.151512 0.615628 0.102768 0.280443 -1.298268 -0.660583 -0.038577
175081 -0.037588 -0.344203 1.648160 0.043454 -0.538057 -1.380873 -1.250300 0.175000 0.753466 -0.461483 0.342089
68019 0.294760 -0.592164 1.361861 -0.695486 0.318040 -1.361088 0.853521 -0.512484 -0.406464 2.187382 -0.180787
118082 1.200928 -0.318027 0.013308 -0.859017 1.212511 -1.757468 0.295551 -0.535305 -0.655033 1.456558 -0.454839
In [30]:
X_train.shape
Out[30]:
(213171, 29)

Графический анализ.¶

Сравним матрицы диаграмм рассеяний для признаков X_train до и после компонентного анализа (PCA), пространство признаков при этом сократилось с 29 до 11 измерений и желательно оценить насколько такое сжатие позволило выделить границы, подчеркнуть зависимости, сократить разброс в данных. Это позволит понять, каким образом значения распределяются по соотв. классам и возможно ли выделить кластеры по тем или иным признакам (такие признаки скорее всего будут наиболее полезны).

In [31]:
# сохраним преобразованную матрицу признаков для графич. анализа
train_pca_plot = pca_united_scaled.fit_transform(X_train)
train_pca_plot.shape
Out[31]:
(213171, 11)
In [32]:
%matplotlib inline

# выведем матрицу диаграмм рассеяния
print('Матрица диаграмм рассеяния набора данных до компонентного анализа' 
      + f' (для {X_train.shape[1]} признаков):')

## обрежем записи до 50 000 штук, чтобы не получить слишком плотное закрашивание
## записи при train_test_split были перемешаны, получится 50 000 случайных примеров из выборки
sns.pairplot(
    pd.concat([X_train, y_train], axis=1)[:50000]
    , diag_kind='kde'
    , hue='TARGET'
    , corner=True
    , plot_kws=dict(marker=".", linewidth=1, alpha=0.3)).fig.set_size_inches(30,30)
plt.show()

Матрица диаграмм рассеяния набора данных до компонентного анализа (для 29 признаков):
In [33]:
# добавим столбец с целевым признаком в набор данных полученный на выходе из PCA

## обрежем записи до 50-100 000 штук, чтобы не получить слишком плотное окрашивание
## (т.к. записи при train_test_split были перемешаны, 
## то на выходе получится 50-100 000 случайных примеров из выборки)

train_pca_plot = pd.concat(
    [pd.DataFrame(
        train_pca_plot[:100000]
        , columns=range(train_pca_plot.shape[1])
        , index=X_train.index[:100000]
    ), y_train[:50000]]
    , axis=1)

# выведем матрицу диаграмм рассеяния
print('Матрица диаграмм рассеяния набора данных после компонентного анализа' 
      + f' (для {train_pca_plot.shape[1]-1} признаков):')

sns.pairplot(
    train_pca_plot
    , diag_kind='kde'
    , hue='TARGET'
    , corner=True
    , plot_kws=dict(marker=".", linewidth=1, alpha=0.3)).fig.set_size_inches(30,30)
plt.show()

Матрица диаграмм рассеяния набора данных после компонентного анализа (для 11 признаков):

Довольно очевидной становится польза от использования метода главных компонент по крайней мере в целях визуализации и в целях сокращения времени на анализ. Хотя было решено использовать для визуализации плоские графики, они всё-равно получаются более информативными.

Распределения сгладились и в частных случаях стали больше похожи на нормальные, кластеризация на некоторых признаках становится отчётливее, скорее всего линейные (и нелинейные) зависимости при их наличии также проявились бы яснее.

Наряду с этим неитерпретируемость самих признаков после компонентного анализа только намекает, о том, что в исходных данных полезная информация содержалась. Проверка полезности таких преобразований приводится в следующем разделе.

Обучение моделей.¶

Перед началом исследования определимся с применяемыми ML-алгоритмами и метриками оценки эффективности моделей.

Для обучения будем пользоваться линейными и нелинейными моделями:

  • *логистической регрессией*,
  • *случайным лесом*,
  • В качесиве базовой оценки примем *константную модель* (DummyClassifier).

В задаче классификации оперировать только метрикой accuracy будет недостаточно, нужно учитывать также изменения precision и recall. С этим вполне справится F1-мера, но поскольку она отражает среднее гармоническое оценок для конкретно выбранного порога чувствительности, то диапазон возможных изменений (для всех порогов) удобнее оценивать по ROC-AUC.

Константная модель.¶

Получим оценку константной/случайной модели по результатам работы DummyClassifier, примем её за базовую. Последовательно передадим DummyClassifier гиперпараметры strategy='uniform' и strategy='constant', чтобы получить случайную модель и константную, которая отвечает "1" независимо от значений на входе.

In [34]:
%%time

# создадим словарь для хранения значений ROC-AUC мер
roc_auc_scores = {}

dummy_random = DummyClassifier(random_state=state, strategy='uniform')
dummy_const = DummyClassifier(random_state=state, strategy='constant', constant=1)

model = [dummy_random, dummy_const]
name = ['CЛУЧАЙНАЯ','КОНСТАНТНАЯ']

# цикл по 2 моделям
for i in range(len(model)):
    clf = Pipeline([           # Pipeline с предобработкой и моделью
        ('features_prepr', pca_united_scaled),
        ('estimator', model[i])
    ])
    
    # получение оценок моделей
    print(f'модель - {str(model[i])[:16]}...{str(model[i])[-20:-1]})')
    print(f'ROC-AUC мера модели ({name[i]}):',
          cross_val_score(clf,  X_train, y_train, 
                          scoring='roc_auc', cv=5).mean(), '\n') # ROC-AUC после кросс-валидации
    
    # добавим результаты в словарь
    roc_auc_scores[f'{name[i]} модель'] = (
        cross_val_score(clf, X_train, y_train, scoring='roc_auc', cv=5).mean().round(4)
    )

модель - DummyClassifier(... strategy='uniform')
ROC-AUC мера модели (CЛУЧАЙНАЯ): 0.5 

модель - DummyClassifier(...strategy='constant')
ROC-AUC мера модели (КОНСТАНТНАЯ): 0.5 

Wall time: 19.8 s

Логистическая регрессия.¶

Получив base-line решение, сравним его с результатами других моделей. *Логистическая регрессия* в этом списке первая. Результы работы модели со станд. настройками будут ужасающе низкими, подберём гиперпараметры с помощью методов оптимизации.

Лучшие гиперараметры найдём случайным поиском по сетке с кросс-валидацией (RandomizedSearchCV), поскольку он позволяет сильно экономить время на вычислениях:

  • для устранения переобучения предполагается в основном варьировать параметром C: Inverse of regularization strength, чтобы достаточно понизить веса признаков на обучающей выборке,

  • подбирать значение лучшего порога не станем.

In [35]:
%%time

model = LogisticRegression(random_state=state)

# Pipeline с предобработкой и моделью
clf = Pipeline([
    ('features_prepr', pca_united_scaled),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__penalty': ['l2'], #['l1', 'l2', 'elasticnet'],
    'estimator__C': [911000, 1000000], #np.linspace(0.00001, 1000000, 10000),
    #'estimator__max_iter': range(100,1000,100),
    'estimator__class_weight': ['balanced'], #['balanced', 'None'],
    'estimator__solver': ['liblinear'] #['saga', 'liblinear']
}

LR_clf = RandomizedSearchCV(
    clf,
    param_distributions=grid,
    random_state=state,
    n_iter=10,
    scoring='roc_auc',
    cv=4, verbose=1, n_jobs=-1
)

LR_clf.fit(X_train, y_train)

# из истории:
#- Лучший ROC-AUC =  0.67438
#- Лучшие параметры: {'solver': 'liblinear', 'penalty': 'l2', 'class_weight': 'balanced', 'C': 911191.1191127993}

# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', LR_clf.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:LR_clf.best_params_.get(i) for i in LR_clf.best_params_})
print('- Лучшая модель:\n\t доступна из "LR_clf.best_estimator_"\n')

# добавим результаты в словарь
roc_auc_scores[f'Логистическая регрессия'] = LR_clf.best_score_

Fitting 4 folds for each of 2 candidates, totalling 8 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.68012
- Лучшие параметры:
	 {'solver': 'liblinear', 'penalty': 'l2', 'class_weight': 'balanced', 'C': 911000}
- Лучшая модель:
	 доступна из "LR_clf.best_estimator_"

Wall time: 8.34 s

Модель, построенная на основе логистической регрессии, дала прирост качества относительно случайного и константного классификаторов, но вполне успешным такой результат признать трудно, поскольку метрики получились плохими.

  • Модель линейная, ей сложно построить сигмоиду в таком многомерном пространстве для убедительного разделения объектов по классам.

  • Наврядли иная обработка признаков дала бы какой-то прирост (допустим, с сохранением большего кол-ва признаков или без их сжатия PCA-алгоритмом), а вот время на обработку увеличилось бы значительно из-за пресловутоо "проклятия размерности". Небольшая проверка этого предположения приведена ниже.

  • Вполне возможно, что выставив порог срабатывания отличным от стандартного = 0.5, можно улучшить результат и даже значительно, но в рамках данной работы этот вопрос не рассматривался из-за временных ограничений. Потребуется построить ROC-AUC-кривую и FPR-TPR-кривую, рассмотреть их формы и сравнить оценки классификатора на разных порогах через метод predict_proba().

In [36]:
%%time

## Альт.вариант логистич.регрессии:
## (РСА на всех признаках, нет итогового масштабирования всех признаков)

model = LogisticRegression(random_state=state)

# новый Pipeline с предобработкой и моделью
clf_1 = Pipeline([
    ('features_prepr_pca_nonsclaed', prepr_union),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__penalty': ['l1', 'l2'], #['l1', 'l2', 'elasticnet'],
    'estimator__C':  [0.001, 1, 10000, 472600], #np.linspace(0.00001, 1000000, 10000),
    #'estimator__max_iter': range(100,1000,100),
    'estimator__class_weight': 'balanced', #['balanced', 'None'],
    'estimator__solver': ['saga', 'liblinear']
}

# новый классификатор
LR_clf_1 = RandomizedSearchCV(
    clf_1,
    param_distributions=grid,
    random_state=state,
    n_iter=20,
    scoring='roc_auc',
    cv=4, verbose=1, n_jobs=-1
)

LR_clf_1.fit(X_train, y_train)

# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', LR_clf_1.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:LR_clf_1.best_params_.get(i) for i in LR_clf_1.best_params_})
print('- Лучшая модель:\n\t доступна из "LR_clf_1.best_estimator_"\n')

# добавим результаты в словарь
#roc_auc_scores[f'Логистическая регрессия'] = LR_clf_1.best_score_

Fitting 4 folds for each of 20 candidates, totalling 80 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.61903
- Лучшие параметры:
	 {'solver': 'saga', 'penalty': 'l1', 'class_weight': 'a', 'C': 10000}
- Лучшая модель:
	 доступна из "LR_clf_1.best_estimator_"

Wall time: 1min 28s
In [37]:
%%time

## Альт.вариант логистич.регрессии:
## (РСА+нормализация на колич. признаках, категориальные только перекодированы)

# новый трансформер для объединения категориальных и числовых признаков после обработки
prepr_union_1 = FeatureUnion([
    ('numeric', num_pca_pipe)
    ,('categorical_nonpca', ct_cat)
])

# новый Pipeline с предобработкой и моделью
clf_2 = Pipeline([
    ('features_prepr_pca_nonsclaed', prepr_union_1),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__penalty': ['l1'],
    'estimator__C':  [123000], #np.linspace(1000, 200000, 200),
    #'estimator__max_iter': range(100,1000,100),
    'estimator__class_weight': ['balanced'],
    'estimator__solver': ['liblinear']
}

LR_clf_2 = RandomizedSearchCV(
    clf_2,
    param_distributions=grid,
    random_state=state,
    n_iter=1,
    scoring='roc_auc',
    cv=4, verbose=1, n_jobs=-1
)

LR_clf_2.fit(X_train, y_train)

# из истории:
#- Лучший ROC-AUC =  0.69631
#- Лучшие параметры: {'solver': 'liblinear', 'penalty': 'l1', 'class_weight': 'balanced', 'C': 122912.2912378938}

#- Лучший ROC-AUC =  0.70005
# Лучшие параметры: {'solver': 'liblinear', 'penalty': 'l2', 'class_weight': 'balanced', 'C': 100}


# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', LR_clf_2.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:LR_clf_2.best_params_.get(i) for i in LR_clf_2.best_params_})
print('- Лучшая модель:\n\t доступна из "LR_clf_2.best_estimator_"\n')

# добавим результаты в словарь
roc_auc_scores[f'Логистическая регрессия_2'] = LR_clf_2.best_score_

Fitting 4 folds for each of 1 candidates, totalling 4 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.69981
- Лучшие параметры:
	 {'solver': 'liblinear', 'penalty': 'l1', 'class_weight': 'balanced', 'C': 123000}
- Лучшая модель:
	 доступна из "LR_clf_2.best_estimator_"

Wall time: 2min 30s

Использование RandomizedSearchCV для настройки гиперпараметров модели вместо перебора всех вариантов по сетке GridSearchCV даёт заметный и прогнозируемый выигрыш в скорости выполнения. Это по-прежнему удивительно, а благодаря теории вероятностей и ЦПТ результаты вдобавок получаются близкими по качеству.

K-ближайших соседей.¶

Обучение модели с помощью алгоритма Взвешенных K-ближайших соседей изначально не рассматривалось по причине величины выборки, но было выполнено "из спортивного интереса" в сокращённом формате.

Основной проблемой является трудоёмкость вычислений для алгоритмов, основанных на метрических расстояниях. Результаты, может, и будут достигаться лучшие, но простое запоминание и расчёт расстояний для каждого объекта в случае большой выборки крайне замедлит процесс получения результатов и будут очень ресурсоёмкими. В качестве выхода можно пробовать дальнейшее снижение размерности, но подробно этот вопрос не исследовался в данной работе.

In [38]:
%%time

model = KNeighborsClassifier()

# новый Pipeline с предобработкой и моделью, признаки дополнительно сжимаются

## повторное сжатие PCA
add_pca_pipe = Pipeline([('pca_united_scaled', pca_united_scaled),
                         ('add_PCA', PCA(0.8))])
# Pipeline
clf = Pipeline([
    ('features_prepr_pca_plus', add_pca_pipe),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__n_neighbors': [40,100, 200], #range(10,41,5),
    'estimator__weights': [lambda x: 1/(x+1)], #['distance', lambda x: 1/(x+1)],
    'estimator__algorithm': ['auto'],
    'estimator__leaf_size': [10, 30], #range(30,101,15),
    'estimator__p': [2]
}

KNС_clf = RandomizedSearchCV(
    clf,
    param_distributions=grid,
    random_state=state,
    n_iter=6,
    scoring='roc_auc',
    cv=5, verbose=1, n_jobs=-1
)

KNС_clf.fit(X_train, y_train)

# из истории:
#- Лучший ROC-AUC =  0.64263
#- Лучшие параметры: {'weights': 'distance', 'p': 2, 'n_neighbors': 50, 'leaf_size': 65, 'algorithm': 'auto'}

# выведем результаты на экран:
print('\nK-ближайших соседей\n')
print('- Лучший ROC-AUC = ', KNС_clf.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:KNС_clf.best_params_.get(i) for i in KNС_clf.best_params_})
print('- Лучшая модель:\n\t доступна из "KNС_clf.best_estimator_"\n')

# добавим результаты в словарь
roc_auc_scores[f'Взвешенный K-ближайших соседей'] = KNС_clf.best_score_

Fitting 5 folds for each of 6 candidates, totalling 30 fits

K-ближайших соседей

- Лучший ROC-AUC =  0.67339
- Лучшие параметры:
	 {'weights': <function <lambda> at 0x00000281C155A8B0>, 'p': 2, 'n_neighbors': 200, 'leaf_size': 10, 'algorithm': 'auto'}
- Лучшая модель:
	 доступна из "KNС_clf.best_estimator_"

Wall time: 5min 14s

Случайный лес¶

Бустинговые и нейросетевые алгоритмы в данной работе не рассматривались и для поиска вариантов улучшения остаётся воспользоваться ансамблевыми методами. В данном случае *случайным лесом*.

Хотя бэггинговый классификатор над деревьями и обладает явными достоинствами в виде нечувствительности к масштабам данных, небольшим числом значимых гиперпараметров и способностью воспроизводить закономерности данных, но ресурсоёмкоесть вычислений на выборке подобного объёма будет велика (по крайней мере для домашнего ПК автора).

  • Вдобавок, из-за PCA-алгоритма и получающейся преобразованной выборки, про интерпретируемость результатов говорить не приходится. Что при хорошем качестве предсказаний можно было бы принять как данность.

  • Дождаться окончания обучения простой модели (со станд. настройками ГП) и на полноразмерной тренировочной выборке за приемлемое время при первом запуске не удалось вовсе, поэтому были предприняты попытки для уменьшения влияния большой размерности.

  • В итоге модели обучались на нескольких выборках для поиска баланса между временем расчёта и качеством:

    • на сокращённой и полной выборке

    • с признаками в исходном виде (без сжатия, без масштабирования, с порядковым кодированием переменных вместо прямого, т.к. деревья прекрасно работают и с ним) так и с обработанной выборкой (сжатие, двойное сжатие, масштабирование).

      как результативный был принят подход обучать модель на выборке pca_united_scaled (т.е. 10 признаков после PCA и масштабирования)

  • Отдельно хочется отметить ценность кросс-валидации при поиске ГП как инструмента, позволяющего получать адекватные и устойчивые оценки качества моделей уже на тренировочной выборке. Без кросс-валидации модели легко переобучаются на всём объёме данных и на тестовой выборке делают предсказания близкие к случайным.

In [39]:
%%time

model = RandomForestClassifier(random_state=state)

## Базовый вариант со случ.лесом:
## (РСА и нормализации нет, колич.признаки в исходном виде, категориальные-перекодированы)

# новый трансформер для категориальных и числовых признаков после обработки
ct_cat_num = ColumnTransformer([
    ("cat_preprocess"
     , OrdinalEncoder(handle_unknown='use_encoded_value', unknown_value = -1)
     , cat_features)
], remainder='passthrough')

# новый Pipeline с предобработкой и моделью
clf_rf = Pipeline([
    ('features_prepr_non_pca_non_scaled', ct_cat_num),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__n_estimators': [300], #range(100, 601, 20),
    'estimator__criterion': ['gini'],
    'estimator__max_depth': [9, 13], #range(5, 41, 2),
    'estimator__class_weight': ['balanced'], #['balanced', None],
    'estimator__max_leaf_nodes': [100] #range(1, 101, 10)
}

RF_clf = RandomizedSearchCV(
    clf_rf,
    param_distributions=grid,
    random_state=state,
    n_iter=1,
    scoring='roc_auc',
    cv=5, verbose=1, n_jobs=-1
) 

RF_clf.fit(X_train, y_train)

# из истории:
#- Лучший ROC-AUC =  0.693
#- Лучшие параметры: {'n_estimators': 500, 'max_depth': 50, 'criterion': 'gini', 'class_weight': 'balanced'}

# выведем результаты на экран:
print('\nСлучайный лес\n')
print('- Лучший ROC-AUC = ', RF_clf.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:RF_clf.best_params_.get(i) for i in RF_clf.best_params_})
print('- Лучшая модель:\n\t доступна из "RF_clf.best_estimator_"\n')

# добавим результаты в словарь
roc_auc_scores[f'Случайный лес'] = RF_clf.best_score_

Fitting 5 folds for each of 1 candidates, totalling 5 fits

Случайный лес

- Лучший ROC-AUC =  0.7009
- Лучшие параметры:
	 {'n_estimators': 300, 'max_leaf_nodes': 100, 'max_depth': 13, 'criterion': 'gini', 'class_weight': 'balanced'}
- Лучшая модель:
	 доступна из "RF_clf.best_estimator_"

Wall time: 2min 52s
In [40]:
%%time

model = RandomForestClassifier(random_state=state)

## Альт.вариант со случайным лесом:
## (РСА и нормализации на всех признаках (однократно, на колич. и категор. по-отдельности)

# первоначальный Pipeline с предобработкой (кодирование+РСА+нормализация) и моделью
clf = Pipeline([
    ('features_prepr', pca_united_scaled),
    ('estimator', model)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid = {
    'estimator__n_estimators': [500], #range(400, 501, 20),
    'estimator__criterion': ['gini'],
    'estimator__max_depth': [7, 13], #range(20, 71, 10),
    'estimator__class_weight': ['balanced'], #['balanced', None],
    'estimator__max_leaf_nodes': [101] #range(1, 101, 2)
    #'ccp_alpha': range(0, 1000, 10)
}

RF_clf_1 = RandomizedSearchCV(
    clf,
    param_distributions=grid,
    random_state=state,
    n_iter=2,
    scoring='roc_auc',
    cv=4, verbose=1, n_jobs=-1
) 

RF_clf_1.fit(X_train, y_train)

#print(accuracy_score(y_train[:100000], clf_rf.predict(X_train[:100000])))
#print(f1_score(y_train[:100000], clf_rf.predict(X_train[:100000])))
#print(roc_auc_score(y_train[:100000], clf_rf.predict(X_train[:100000])))

# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', RF_clf_1.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:RF_clf_1.best_params_.get(i) for i in RF_clf_1.best_params_})
print('- Лучшая модель:\n\t доступна из "RF_clf.best_estimator_"\n')

# добавим результаты в словарь
roc_auc_scores[f'Случайный лес_2'] = RF_clf_1.best_score_

Fitting 4 folds for each of 2 candidates, totalling 8 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.68506
- Лучшие параметры:
	 {'n_estimators': 500, 'max_leaf_nodes': 101, 'max_depth': 13, 'criterion': 'gini', 'class_weight': 'balanced'}
- Лучшая модель:
	 доступна из "RF_clf.best_estimator_"

Wall time: 8min 9s

  • Результаты иногда становятся лучше достижений логистической регрессии, но их получение было значительно более трудоёмкое. Хотя результат и получился довольно посредственный, остановимся тем не менее на имеющихся, чтобы по крайней мере можно было завершить исследование.

  • Модель явно можно улучшить подбором гиперпараметров, а также настройкой порога классификации, но эта деятельность в рамки учебной работы не вошла из-за временных и ресурсных ограничений.

Результаты на тестовой выборке, выводы, альтернативы.¶

Сравним оценки моделей на отложенной тестовой выборке и составим выводы с использованием матрицы ошибок/отчётов о классификации и значимости признаков.

In [41]:
# соберём результаты в общий датафрейм

# ROC-AUC построенных моделей
print('ROC-AUC построенных моделей:')
roc_auc_df = (
    pd.DataFrame.from_dict(roc_auc_scores,
                           orient='index',
                           columns=['ROC-AUC_train+cv'])
#    .sort_values(by='ROC-AUC_train+cv', ascending=True)
    .round(4)
)

# соберём ROC-AUC моделей на тестовой выборке и добавим его столбцом в датафрейм roc_auc_df
dummy = (
    Pipeline([
        ('prepr', pca_united_scaled),
        ('estimator', DummyClassifier(random_state=state, strategy='uniform'))
    ])
    .fit(X_train, y_train)
)

roc_auc_test = (                                                        ##ROC-AUC-score
    [roc_auc_score(y_test, dummy.predict(X_test))                       ## случайной модели
     , 0.5                                                              ## константной модели
     , roc_auc_score(y_test, KNС_clf.best_estimator_.predict(X_test))   ## K-ближайших
     , roc_auc_score(y_test, LR_clf.best_estimator_.predict(X_test))    ## логистич.регрессии
     , roc_auc_score(y_test, LR_clf_2.best_estimator_.predict(X_test))  ## логистич.регрессии_2
     , roc_auc_score(y_test, RF_clf.best_estimator_.predict(X_test))    ## cлучайного леса
     , roc_auc_score(y_test, RF_clf_1.best_estimator_.predict(X_test))] ## cлучайного леса_2
)

roc_auc_df['ROC-AUC_test'] =  pd.Series(roc_auc_test, index=roc_auc_df.index)
roc_auc_df.sort_values(by='ROC-AUC_train+cv', ascending=True)

ROC-AUC построенных моделей:
Out[41]:
ROC-AUC_train+cv ROC-AUC_test
CЛУЧАЙНАЯ модель 0.5000 0.498070
КОНСТАНТНАЯ модель 0.5000 0.500000
Взвешенный K-ближайших соседей 0.6734 0.642263
Логистическая регрессия 0.6801 0.500000
Случайный лес_2 0.6851 0.628132
Логистическая регрессия_2 0.6998 0.622773
Случайный лес 0.7009 0.641955

Выводы:

  • Оценки на тестовой выборке снилизись, что естественно.

  • Видно, что случайный лес и кластеризацию К-ближайших тем или иным способом (довольно случайно) удалось обучить лучше остальных, но доверять этому результату не хочется, из-за времени работы алгоритмов. Особенно сложным в использовании кажется кластерный алгоритм, поскольку он по сути не имеет стадии "обучения" для создания обобщающих правил и по своей природе просто должен хранить в памяти целиком всю выборку для рассчёта расстояний.

  • В качестве лучшего на данном этапе признаем результат логистической регрессии, он кажется наиболее устойчивым.

  • Большой объём данных сильно сказывается на скорости обучения моделей и добиться хорошего результата за счёт отбрасывания признаков и уменьшения размерности не слишком просто, поскольку выигрыш в скорости обучения компенсируется затратами времени на соотв. подготовку данных (препроцессинг и трансформирование).

  • Настроить алгоритмы кажется вполне реальной задачей, изменяя порог срабатывания моделей через predict_proba(), но в рамках данной работы этот приём не применялся, целью было изучение методов снижения размерности и практика использования конвейеров и трансформеров.

  • Возникает подозрение, что при обилии признаков и объектов и избыточности в информации проще сразу использовать бустинговые модели, которые обучаются с компенсацией ошибок или нейронные сети (правда насколько мне известно для табличных данных каких-то общепринятых архитектур не существует), чем подбирать ГП, т.к. этот процесс трудоёмок и требует времени.

Значимость признаков и модификация лучшей модели.¶

В данной работе нет принципиального требования улучшить результаты, поэтому просто дополнительно взглянем на значимость признаков полученной модели и проанализируем матрицу ошибок классификации на тестовой выборке.

In [42]:
# построим график значимости признаков на основе permutation_importance
features_names = X_test.columns #названия столбцов

# получение словаря с permutation_importance
result = (
    permutation_importance(
        LR_clf_2.best_estimator_
        , X_test, y_test
        , n_repeats=30
        , random_state=state)
)

# преобразуем в pd.Series для построения гистограммы
model_importances = pd.Series(result['importances_mean'], index=features_names)

# построение гистограммы значений со станд. отклонением
fig, ax = plt.subplots(figsize=(10,10))
model_importances.plot.bar(yerr=result['importances_std'], ax=ax, rot=90)
ax.set_title("Feature importances")
ax.set_ylabel("Среднее значение")
ax.set_xlabel("Features")
fig.tight_layout()
plt.show()

Вывод:

  • Оценка качества моделей в принципе невысока, выдающихся результатов достигнуто не было, значимость признаков неравномерна и не слишком велика.

  • можно видеть, что далеко не все признаки оказались полезными для предсказаний модели. Некоторые и вовсе вносят отрицательный вклад в результат. В дальнейшем можно пробовать упростить модель, исключив малозначимые признаки и признаки с отрицательным вкладом (либо попробовать учитывать их с обратным знаком).

  • Но есть предположение, что наличие категориальных признаков и использование PCA, преобразующего значения нелинейно, обойтись изменением знака, описанными выше, не позволят. Хотя широкое использование FunctionTransformer и make_column_selector для создания своих конвейеров и трансформеров обработки признаков меня очень впечатлило получающейся в результате гибкостью.

In [43]:
# удалим отрицательно- и малозначимые признаки
# подготовим список соотв. столбцов
mi_minus = model_importances[model_importances<=0.001].index
mi_minus
Out[43]:
Index(['FLAG_OWN_REALTY', 'CNT_CHILDREN', 'AMT_CREDIT', 'NAME_TYPE_SUITE',
       'NAME_INCOME_TYPE', 'NAME_EDUCATION_TYPE', 'NAME_HOUSING_TYPE',
       'DAYS_BIRTH', 'DAYS_REGISTRATION', 'DAYS_ID_PUBLISH', 'FLAG_WORK_PHONE',
       'FLAG_CONT_MOBILE', 'FLAG_PHONE', 'FLAG_EMAIL', 'REGION_RATING_CLIENT',
       'HOUR_APPR_PROCESS_START', 'REG_REGION_NOT_LIVE_REGION',
       'REG_CITY_NOT_WORK_CITY', 'ORGANIZATION_TYPE',
       'OBS_30_CNT_SOCIAL_CIRCLE', 'DEF_30_CNT_SOCIAL_CIRCLE',
       'DAYS_LAST_PHONE_CHANGE', 'FLAG_DOCUMENT_3'],
      dtype='object')
In [44]:
# признаки удалим вручную
X_train_mod = X_train.drop(columns=mi_minus).copy()
X_test_mod = X_test.drop(columns=mi_minus).copy()
In [45]:
%%time

## Альт.вариант логистич.регрессии:
## (сокращённая после получения оценок значимости признаков)

model = LogisticRegression(random_state=state)

# оптимизация гиперпарметров по сетке "родительской" модели
grid = {
    'estimator__penalty': ['l1'],
    'estimator__C': np.linspace(104000, 125000, 20),
    #'estimator__max_iter': range(100,1000,100),
    'estimator__class_weight': ['balanced'],
    'estimator__solver': ['liblinear']
}

LR_clf_2_mod = RandomizedSearchCV(
    clf_2,
    param_distributions=grid,
    random_state=state,
    n_iter=10,
    scoring='roc_auc',
    cv=4, verbose=1, n_jobs=-1
)

LR_clf_2_mod.fit(X_train_mod, y_train)

# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', LR_clf_2_mod.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:LR_clf_2_mod.best_params_.get(i) for i in LR_clf_2_mod.best_params_})
print('- Лучшая модель:\n\t доступна из "LR_clf_2_mod.best_estimator_"\n')

Fitting 4 folds for each of 10 candidates, totalling 40 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.67151
- Лучшие параметры:
	 {'solver': 'liblinear', 'penalty': 'l1', 'class_weight': 'balanced', 'C': 105105.26315789473}
- Лучшая модель:
	 доступна из "LR_clf_2_mod.best_estimator_"

Wall time: 26.3 s
In [46]:
# оценка модифицированной "облегчённой" модели на тесте
roc_auc_score(y_test, LR_clf_2_mod.predict(X_test_mod))
Out[46]:
0.6192376678727323

Несколько жертвуя качеством можно продолжить сокращать кол-во используемых признаков, скорость работы моделей при этом возрастёт за счёт снижения качества.

Матрица ошибок лучшей модели.¶

Ниже взглянем на матрицу ошибок лучшей модели.

In [47]:
# classification_report
y_true = y_test
y_pred = LR_clf_2.best_estimator_.predict(X_test)
target_names = ['class 0', 'class 1']
print('Матрица мер классов целевого признака (на тестовой выборке):\n')
print(classification_report(y_true, y_pred, target_names=target_names))
#print(confusion_matrix(y_true, y_pred))

# матрица ошибок
print('\nМатрица ошибок (на тестовой выборке):')
ConfusionMatrixDisplay(confusion_matrix(y_true, y_pred), display_labels=target_names).plot()
plt.show()

Матрица мер классов целевого признака (на тестовой выборке):

              precision    recall  f1-score   support

     class 0       0.95      0.66      0.78     83960
     class 1       0.14      0.62      0.23      7400

    accuracy                           0.66     91360
   macro avg       0.55      0.64      0.50     91360
weighted avg       0.89      0.66      0.74     91360


Матрица ошибок (на тестовой выборке):

В представленной модели результаты не видятся слишком удачными. Выдающееся количество ложно-положительных срабатываний, явно говорит о путях улучшения классификации в дальнейшем в зависимости от ценности тех или иных клиентов для бизнеса. Здесь просматривается возможность изменения порога срабатывания классификации (thresholod) для варьирования метриками precision и recall. При наличии свободных ресурсов эти изменения могут быть реализованы в доп. разделе.

Вывод:

  • Дисбаланс классов влияет на результаты предсказаний модели. Разумным копмромиссом на фоне снижения размерности признакогово пространства было бы масштабирование объектов целевых классов для снижения влияния дисбаланса, например увеличение кол-ва объектов младшего класса или взвешивание предсказаний моделей.

  • Зная о преимущественном характере ложных срабатываний модели, можно дополнительно скорректировать её поведение различными способами (сдвига порога, взвешивание классов при обучении, подбор ГП). В данной работе эти приёмы не затронуты в полной мере в силу демонстрационного характера для других задач.

  • При подготовке признаков сокращение размерности и нормализация числовых значений позволяют хотя бы приступить к этапу обучения моделей на локальной машине, что уже можно признать полезным. Большие наборы данных имеет смысл подготавливать именно способом уменьшения размерности признакового пространства, а также скорее всего использовать бустинговые или даже перцептронные модели.

  • По всей видимости для НЕ-нейронных моделей будет всегда иметься компромисс между выбором большего ROC_AUC и между статистиками классфикации "precision-recall" и выбор стоит делать опираясь на требования бизнеса.

Некоторые альтернативы и компромиссы.¶

После получения оценок на тестовой выборке можно провести обучение лучшей модели на объединённой выборке, но из-за низкого итогового качества этот пункт, пожалуй, можно опустить. Ниже бегло взглянем на альтернативный подход к решению с использованием других ML-моделей и остановимся в изысканиях.

Бустинг.¶

In [48]:
# пробуем бустинг
from catboost import CatBoostClassifier
from catboost import cv, Pool
In [49]:
## пробный запуск бустинга
## используем классификатор CatBoost из-за его возможностей по работе с категориальными переменными

# инициализация модели, устанавливаем ГП
model_cat = CatBoostClassifier(
    loss_function='Logloss'
    , eval_metric='F1'
    , iterations=500
    , custom_loss=['AUC', 'Accuracy']
    , random_seed = 42
    , learning_rate = 2
    , bagging_temperature=1
    , random_strength=10
    #, thread_count=3
    , l2_leaf_reg = 100
    #, save_snapshot=True
    #, snapshot_file='snapshot_best.bkp'
    #, od_type='Iter'
    #, od_wait=20
    , use_best_model=True
)

# обучение на тренировочной выборке с валидацией по тестовой чтобы найти лучшую модель
model_cat.fit(
    X_train, y_train
    , cat_features=[col for col in X_train.columns
                    if X_train[col].dtype != 'float64' and len(X_train[col].unique())<60]
    , eval_set=(X_test, y_test)   # eval_set, чтобы использовать use_best_model
    #, logging_level='Silent'
    , plot=True
    , verbose=50
    , early_stopping_rounds=20
    , use_best_model=True
)

MetricVisualizer(layout=Layout(align_self='stretch', height='500px'))
learning rate is greater than 1. You probably need to decrease learning rate.
learning rate is greater than 1. You probably need to decrease learning rate.
learning rate is greater than 1. You probably need to decrease learning rate.

0:	learn: 0.0000000	test: 0.0000000	best: 0.0000000 (0)	total: 464ms	remaining: 3m 51s
50:	learn: 0.0003472	test: 0.0005402	best: 0.2274433 (35)	total: 10.7s	remaining: 1m 34s
Stopped by overfitting detector  (20 iterations wait)

bestTest = 0.2274433023
bestIteration = 35

Shrink model to first 36 iterations.
Out[49]:
<catboost.core.CatBoostClassifier at 0x281c01e9ca0>
In [50]:
# roc_auc_score на тестовой выборке
roc_auc_score(y_test, model_cat.predict(X_test))
Out[50]:
0.6030520324993882
In [51]:
# f1_score на тестовой выборке
f1_score(y_test, model_cat.predict(X_test))
Out[51]:
0.22744330233317014

Для получения хоть какого-то результата на выборке уменьшенной размерности (теже X_train и X_test на 30 признаков, что и у моделей из предыдущего раздела) в модель пришлось передавать весьма "агрессивные" значения гиперпараметров. По всей видимости для CatBoost эти признаки оказались далеко не такими информативными как для классич. алгоритмов.

Проверим, улучшится ли качество бустинговой модели, при передаче в неё признаков без изъятий (кроме тех, где было много пропусков). Делать проверку на исходных данных в полном объёме не станем (для экономии времени). Бэггинг над бустингами также не проводится.

In [52]:
# новое разбиение на выборки, пространство признаков без сокращений
X_ctb_train, X_ctb_test, y_ctb_train, y_ctb_test = (
    train_test_split(df.drop(columns=['SK_ID_CURR', 'TARGET']), df['TARGET'], train_size = 0.7
                     , random_state=state, shuffle=True
                     , stratify=df['TARGET'])
)

X_ctb_train.shape, X_ctb_test.shape, y_ctb_train.shape, y_ctb_test.shape
Out[52]:
((213171, 52), (91360, 52), (213171,), (91360,))
In [53]:
## запуск бустинга для расшир. пространства признаков с новыми ГП

# инициализация модели, устанавливаем ГП
model_cat_full = CatBoostClassifier(
   loss_function='Logloss'
    , eval_metric='F1'
    , iterations=500
    , custom_loss=['AUC', 'Accuracy']
    , random_seed = 42
    , learning_rate = 2
    , bagging_temperature=1
    , random_strength=10
    #, thread_count=3
    , l2_leaf_reg = 100
    #, save_snapshot=True
    #, snapshot_file='snapshot_best.bkp'
    #, od_type='Iter'
    #, od_wait=20
    , use_best_model=True
)

# обучение на тренировочной выборке с валидацией по тестовой чтобы найти лучшую модель
model_cat_full.fit(
    X_ctb_train, y_ctb_train
    , cat_features=[col for col in X_ctb_train.columns
                    if X_ctb_train[col].dtype != 'float64' and len(X_ctb_train[col].unique())<60]
    , eval_set=(X_ctb_test, y_ctb_test)
    #, logging_level='Silent'
    , plot=True
    , verbose=50
    , early_stopping_rounds=50
    , use_best_model=True
)

learning rate is greater than 1. You probably need to decrease learning rate.

MetricVisualizer(layout=Layout(align_self='stretch', height='500px'))
learning rate is greater than 1. You probably need to decrease learning rate.
learning rate is greater than 1. You probably need to decrease learning rate.

0:	learn: 0.0000000	test: 0.0000000	best: 0.0000000 (0)	total: 290ms	remaining: 2m 24s
50:	learn: 0.0042291	test: 0.0018774	best: 0.2252610 (17)	total: 12.5s	remaining: 1m 49s
100:	learn: 0.0045709	test: 0.0016088	best: 0.2268068 (51)	total: 25.4s	remaining: 1m 40s
Stopped by overfitting detector  (50 iterations wait)

bestTest = 0.2268068447
bestIteration = 51

Shrink model to first 52 iterations.
Out[53]:
<catboost.core.CatBoostClassifier at 0x281bfcb5c70>
In [54]:
# roc_auc_score на тестовой выборке
roc_auc_score(y_ctb_test, model_cat_full.predict(X_ctb_test))
Out[54]:
0.5982303992892367
In [55]:
# f1_score на тестовой выборке
f1_score(y_ctb_test, model_cat_full.predict(X_ctb_test))
Out[55]:
0.22680684467514184

Изменения малозначительные, можно остнавиться на предыд. этапе. Судя по всему, бустинг требует отличной от предложенной ранее предобработки и другого подхода при отборе признаков.

В данном исследовании остановим поиск решений среди бустинговых алгоритмов в виду ограниченности времени на исследование, с развитием данного проекта пункт о поиске более удачного решения среди бустингов будет одним из первых в очереди на совершенствование.

## Код для получения "скользящей" оценки среднего качества бустинга (не запускался)

model = ...

params = {
    'loss_function': 'AUC',
    'iterations': 500,
    'custom_loss': 'AUC',
    'random_seed': 42,
    'learning_rate': 0.05
}

cv_data = cv(
    params=params,
    pool=Pool(X_train, label=y_train, cat_features=X_train.select_dtypes(exclude='float64').columns.values),
    fold_count=5, # Разбивка выборки на 5 кусочков
    shuffle=True, # Перемешаем наши данные
    partition_random_seed=0,
    plot=True,    # Визуализация
    stratified=True, 
    verbose=100,
    early_stopping_rounds=20
)

Доп. логистическая регрессия.¶

Исследуем также работоспособность логистической регрессии при более тщательном отборе признаков c добавлением прямого кодирования (OHE) вместо LabelEncoding, поскольку в задачах классификации с потенциально линейно-разделимыми признаками эта модель при должной настройке должна показывать хорошие результаты при высокой производительности.

К сожалению, из-за кратного роста числа признаков после OHE-кодирования, на обучение моделей будет затрачиваться на порядок больше времени, поэтому проверим предположение только в применении к томуже набору признаков X_train, что подавался для обучение моделей выше. На болшем количестве сырых "данных" (вроде X_ctb_train, полученных до сжатия и отбрасывания) время обработки вырастет ещё больше и вдобавок навряд ли даст прирост метрик качества (хоть и в ограниченном объёме, но эта версия всё-таки подвергалась проверке).

In [56]:
%%time
## итоговая попытка улулчшить результат работы моделей (с более тщательным подходом к отбору признаков)
## категориальными признаем не все признаки, отн. к типу int64 и object, а только отвечающие условиям

## условия для выбора категориальных и числовых данных
cat_features_lrf = [col for col in X.columns if X[col].dtype != 'float64' and len(X[col].unique())<60]
num_features_lrf = [col for col in X.columns if not (X[col].dtype != 'float64' and len(X[col].unique())<60)]

## трансформеры для категориальных и числовых данных (параллельные)

# стандартизация числовых признаков
ct_num = ColumnTransformer([
    ("num_preprocess", StandardScaler(), num_features_lrf)
], remainder='drop')

# преобразование данных из соотв. трансформмеров и сжатие PCA
num_pca_pipe = Pipeline([('num_prepr', ct_num), ('num_PCA', PCA(0.95))])

# кодирование категориальных признаков
ct_cat = ColumnTransformer([
    ("cat_preprocess"
     , OneHotEncoder(drop='if_binary', handle_unknown='ignore', sparse=True)
     , cat_features_lrf)
], remainder='drop')


## Псоледний альт.вариант логистич.регрессии:
## (РСА+нормализация на колич. признаках, категориальные только перекодированы с пом. OHE)

# новый трансформер для объединения категориальных и числовых признаков после обработки
prepr_union_lrf = FeatureUnion([
    ('numerical_scaled_pca', num_pca_pipe)
    ,('categorical_ohe', ct_cat)
])

# инициализация модели
model_alt = LogisticRegression(random_state=state)

# новый Pipeline с предобработкой и моделью
clf_lrf = Pipeline([
    ('features_prepr', prepr_union_lrf),
    ('estimator', model_alt)
])

# поиск лучших гиперпарметров по сетке с кросс-валидацией
grid_lrf = {
    'estimator__penalty': ['elasticnet'], #['l1', 'l2'],
    'estimator__C': [5],
                    #[0.001
                    #, *np.linspace(0.01, 1, 3)
                    #, 50
                    #, *np.linspace(100, 1000, 3)
                    #, 5000
                    #, *np.linspace(10000, 200000, 10)
                    #, 500000],
    'estimator__max_iter': [10000],
    'estimator__class_weight': ['balanced'],
    'estimator__solver': ['saga'],
    'estimator__l1_ratio': [0] #[0, 1, 0.5]
}

LR_clf_f = RandomizedSearchCV(
    clf_lrf,
    param_distributions=grid_lrf,
    random_state=state,
    n_iter=1,
    scoring='roc_auc',
    cv=4, verbose=2, n_jobs=-2
)

LR_clf_f.fit(X_train, y_train)

# из истории:
#- Лучший ROC-AUC =  0.69767
#- Лучшие параметры: {'solver': 'saga', 'penalty': 'elasticnet', 'max_iter': 10000, 'l1_ratio': 0, 'class_weight': 'balanced', 'C': 0.01/0.02}

# выведем результаты на экран:
print('\nЛогистическая регрессия\n')
print('- Лучший ROC-AUC = ', LR_clf_f.best_score_.round(5))
print('- Лучшие параметры:\n\t', {str(i)[11:]:LR_clf_f.best_params_.get(i) for i in LR_clf_f.best_params_})
print('- Лучшая модель:\n\t доступна из "LR_clf_f.best_estimator_"\n')

# добавим результаты в словарь
#roc_auc_scores[f'Логистическая регрессия_2_lrf'] = LR_clf_f.best_score_

Fitting 4 folds for each of 1 candidates, totalling 4 fits

Логистическая регрессия

- Лучший ROC-AUC =  0.69887
- Лучшие параметры:
	 {'solver': 'saga', 'penalty': 'elasticnet', 'max_iter': 10000, 'l1_ratio': 0, 'class_weight': 'balanced', 'C': 5}
- Лучшая модель:
	 доступна из "LR_clf_f.best_estimator_"

Wall time: 20min 29s
In [57]:
# roc_auc + F1 настраиваемой модели и моделей предыд. поколения на тестовой выборке
print('Метрики моделей на тестовой выборке:'
      +'\n\n\t* алт. модели LR_clf_f (логистич. регрессия)'
      +'\n\t\t- ROC-AUC '
      +str(np.round(roc_auc_score(y_test, LR_clf_f.best_estimator_.predict(X_test)),4))
      +'\n\t\t- F1 '
      +str(np.round(f1_score(y_test, LR_clf_f.best_estimator_.predict(X_test)),4)))

print('\n\t* станд. модели LR_clf_2 (логистич. регрессия)'
      +'\n\t\t- ROC-AUC '
      +str(np.round(roc_auc_score(y_test, LR_clf_2.best_estimator_.predict(X_test)),4))
      +'\n\t\t- F1 '
      +str(np.round(f1_score(y_test, LR_clf_2.best_estimator_.predict(X_test)),4)))

print('\n\t* алт. модели RF_clf (случайный лес)'
      +'\n\t\t- ROC-AUC '
      +str(np.round(roc_auc_score(y_test, RF_clf.best_estimator_.predict(X_test)),4))
      +'\n\t\t- F1 '
      +str(np.round(f1_score(y_test, RF_clf.best_estimator_.predict(X_test)),4)))

Метрики моделей на тестовой выборке:

	* алт. модели LR_clf_f (логистич. регрессия)
		- ROC-AUC 0.6439
		- F1 0.2292

	* станд. модели LR_clf_2 (логистич. регрессия)
		- ROC-AUC 0.6423
		- F1 0.2281

	* алт. модели RF_clf (случайный лес)
		- ROC-AUC 0.642
		- F1 0.2318

Уточнения, внесённые в процесс отбора признаков и более тщательный подбор гиперпараметров помогли немного нарастить качество логистической регрессии, поскольку предпочтение ещё на этапе сравнения моделей было отдано этой модели. Используем модифицированную версию как итоговую и проверим на ней способ оптимизации в виде настройки чувствительности.

Выбор порога threshold полученной модели.¶

Известно, что точность предсказаний логистической регрессии может сильно меняться в зависимости от порога чувствительности. Выбрав лучшую модель после сравнении с другими альтернативами, рассмотрим также её поведение при сдвиге порога срабатывания threshold от стандартного значения 0.5.

Оценим оптимальное значение threshold для полученной модели путём максимизации F1-меры в качестве примера (с целью максимально увеличить число найденных True Positive значений). В реальности скорее стоит подбирать порог, исходя из из соображений максимизаци метрик precision или recall, уменьшая количество ложных срабатываний, и в зависимости от их ценности в задаче. Подбор значения произведём следующим образом:

  • последовательно будем выбирать значение порогов из диапазона 0..1,
  • по списку вероятностей положительного класса для соотв. значения порога будем передавать подошедшие по условию значения вероятностей и сформируем новый список целевого признака,
  • передадим новый список предсказаний в функцию f1_score для рассчёта F1-меры,
  • сравним результаты,
  • отобразим результат на графике зависимости F1-score от значений threshold (по значениями из списков f1_list и thrsh_list, полученных в цикле),
  • для большей точности можно получать оценки в кросс-валидации, но в силу размеров датасета от этой идеи можно отказаться.
In [58]:
# найдём значения f1-меры для разных значений порога threshold (без кросс-валидации):

def thrsh(model, description):
    print(f'\nУточнение порога чувствительности модели: {description}.\n')

    best_F1 = 0
    best_threshold = 0
    f1_list = []
    thrsh_list = []
    
    ## подготовим список вероятностей положительного класса для модели
    probabilities_one = model.best_estimator_.predict_proba(X_train)[:,1]

    for threshold in np.arange(0, 1, 0.01):
        ## получим список значений True/False для выбранного threshold, 
        ## сравнив с ним вероятности положительного класса
        predicted = probabilities_one > threshold
        F1_score = f1_score(y_train, predicted) #№ рассчитаем f1-меру, проверив условие
        thrsh_list.append(threshold)
        f1_list.append(F1_score)
    
        if F1_score > best_F1:
            best_F1, best_threshold = F1_score, threshold

    f1 = f1_score(y_train, model.best_estimator_.predict(X_train))
    print('Определёны: best_F1_score =' 
          + str(best_F1.round(4))
          + '| best_threshold ='
          + str(best_threshold.round(4))
         )
    txt = f'F1-score по умолчанию (при threshold - 0.5) = {f1.round(4)}'
    print(txt)
    print('-'*len(txt))

    ## строим график по полученным выше значениям threshold и F1 
    ## из списков thrsh_list, f1_list:
    
    plt.figure(figsize=(16,6))
    plt.plot(thrsh_list, f1_list)

    plt.plot([best_threshold,best_threshold],[0,best_F1+0.05],
             linestyle='dotted',
             linewidth=2) ## вертик.линия с лучшим порогом
    
    plt.plot([0.5, 0.5],[0,best_F1+0.05],
             linestyle='dotted',
             alpha=0.7) ## вертик.линия со станд. порогом

    plt.xlabel('Threshold', fontsize=12)
    plt.ylabel('F1-score', fontsize=12)

    plt.title(f'Зависимость F1-меры от значений порога чувствительности модели',
              fontsize=12)
    plt.legend(['F1-мера',
                f'Лучшее значение порога={best_threshold.round(4)}'
                + f'(F1={best_F1.round(4)})',
                f'Станд.порог=0.5 (F1={f1.round(4)})'],
               loc='lower right')
    plt.show()
    
    ## функция будет возвращать значение лучшего порога чувствительности
    return best_threshold
    
best_threshold = thrsh(LR_clf_f, 'Логистическая регрессия "LR_clf_f"')

Уточнение порога чувствительности модели: Логистическая регрессия "LR_clf_f".

Определёны: best_F1_score =0.2537| best_threshold =0.63
F1-score по умолчанию (при threshold - 0.5) = 0.2322
----------------------------------------------------

Порог из соображений максимизации F1-меры был найден, взглянем как он влияет на значения в матрице ошибок классификации, сравнив её с результатом при неизменённом пороге.

In [59]:
# classification_report
y_true = y_test
y_pred = LR_clf_f.best_estimator_.predict(X_test)
y_pred_best_thrsh = LR_clf_f.best_estimator_.predict_proba(X_test)[:,1] > best_threshold

# матрица ошибок
print('\nМатрица ошибок классификации (на тестовой выборке) с порогом 0,5:')
ConfusionMatrixDisplay(confusion_matrix(y_true, y_pred)
                       , display_labels=target_names).plot()
print('\nМатрица ошибок классификации (на тестовой выборке) с порогом best_threshold:')
ConfusionMatrixDisplay(confusion_matrix(y_true, y_pred_best_thrsh)
                       , display_labels=target_names).plot()
plt.show()

Матрица ошибок классификации (на тестовой выборке) с порогом 0,5:

Матрица ошибок классификации (на тестовой выборке) с порогом best_threshold:

Нахождение значения threshold путём максимизации F1-меры привело к увеличению числа найденных True Positive значений.

Скорее всего в реальности полезнее подбирать порог, исходя из из соображений максимизаци метрик *precision* или *recall*, уменьшая таким образом количество ложных срабатываний по классам 1 или 0 в зависимости от их ценности для решения задачи. Решение о чём стоит принимать, исходя из стратегии заказчика (больший потенциал прибыли при большем риске или меньшая, но более стабильная прибыль при меньшем риске в данном случае).

Напоследок оценим, как изменилась значимость признаков для модели логистической регрессии после модификации.

In [60]:
# построим график значимости признаков на основе permutation_importance
features_names = X_test.columns #названия столбцов

# получение словаря с permutation_importance
result = (
    permutation_importance(
        LR_clf_f.best_estimator_
        , X_test, y_test
        , n_repeats=30
        , random_state=state)
)

# преобразуем в pd.Series для построения гистограммы
model_importances = pd.Series(result['importances_mean'], index=features_names)

# построение гистограммы значений со станд. отклонением
fig, ax = plt.subplots(figsize=(10,10))
model_importances.plot.bar(yerr=result['importances_std'], ax=ax, rot=90)
ax.set_title("Feature importances")
ax.set_ylabel("Среднее значение")
ax.set_xlabel("Features")
fig.tight_layout()
plt.show()

Вывод:

  • Регуляризация весов, увеличение кол-ва итераций для сходимости решения и отбор признаков повлияли на изменение значимости признаков для модели, а также на выравнивание их влияния внутри групп.

  • Учёт этих результатов представляет интерес для дальнейшей доработки модели и её облегчения с потенциалом к улучшению. Эти работы могут быть продолжены при дальнейшем развитии данного проекта наряду с поиском иных моделей для достижения лучших результатов.

  • На данном этапе работу можно считать законченной в достаточной для демонстрационных целей мере: отбор и сокращение признакового пространства произведены по запланированному сценарию, модели обучены и сравнены между собой с поиском лучши гиперпараметров и анализом результов и путей улучшения, выявлены направления развития проекта.